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A HIGH ORDER FINITE DIFFERENCE SCHEME WITH SHARP SHOCK 
RESOLUTION FOR THE EULER EQUATIONS 

MARGOT GERRITSEN AND PELLE OLSSON 


Abstract. We derive a high-order finite difference scheme for the Euler equations that sat- 
isfies a semi-discrete energy estimate, and present an efficient strategy for the treatment of 
discontinuities that leads to sharp shock resolution. The formulation of the semi-discrete 
energy estimate is based on a symmetrization of the Euler equations that preserves the ho- 
mogeneity of the flux vector, a canonical splitting of the flux derivative vector, and the use 
of difference operators that satisfy a discrete analogue to the integration by parts procedure 
used in the continuous energy estimate. Around discontinuities or sharp gradients, refined 
grids are created on which the discrete equations are solved after adding a newly constructed 
artificial viscosity. The positioning of the sub-grids and computation of the viscosity are aided 
by a detection algorithm which is based on a multi-scale wavelet analysis of the pressure 
grid function. The wavelet theory provides easy to implement mathematical criteria to detect 
discontinuities, sharp gradients and spurious oscillations quickly and efficiently. 


1. Introduction 

High-order finite difference schemes are attractive for computing Euler flows that have large smooth 
flow regions. Their simplicity makes coding easy while their high order of accuracy allows for the use 
of coarse grids in the regions where the flow is slowly varying. In multi-dimensional applications this 
reduces both computational costs and memory requirements significantly. In the past, two problems 
prevented wide spread use of these schemes. Firstly, it was not known how to construct boundary dif- 
ference operators of sufficient accuracy that lead to (strictly) stable numerical schemes. Secondly, the 
schemes introduce spurious oscillations at the shock that need to be damped by filtering or addition of 
artificial viscosity. Sharp shock resolution using previous techniques has required considerable tuning. 
In this paper we address both issues. We devise a high-order finite difference scheme that satisfies 
an semi-discrete energy estimate and present an efficient strategy for the treatment of discontinuities 
that leads to sharp shock resolution. 

Stability of initial-boundary value problems can be ensured by the energy method. In general 
the method consists of an integration-by-parts procedure which introduces boundary terms that 
are subsequently bounded or eliminated to obtain the desired energy estimate. In [9, 17], Kreiss, 
Scherer and Olsson developed a semi-discrete analogue of this procedure for linear equations and 
special finite difference operators. Olsson showed that extension of the energy method to nonlinear 
systems is possible [18]. He employs a splitting of the flux derivative vector after symmetrization of 
the system. To facilitate the derivation of the energy estimates for the Euler equations, the variable 
transformation should preserve the homogeneity of the flux vector. We show that these conditions 
can be met for a specific class of transformations. 
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Inspired by the stability results we devise a high-order scheme using the aforementioned trans- 
formation and splitting of the flux derivative vector. It results in a system of ordinary differential 
equations that is solved using a TVD Runge Kutta time integration. We refer to it as the Split 
High-Order Entropy-Conserving (SHOEC) scheme. 

In the second and main part of the paper we present an efficient strategy for the treatment of 
discontinuities. Our purpose is twofold: 

• Control spurious oscillations without sacrificing the overall 
computational efficiency of high-order methods. 

• Resolve shocks sharply without having to adjust parameters. 

In achieving these goals the stability properties of the method should not be destroyed. Our strategy 
is based on local grid adaptation. The smooth parts of the flow are resolved on coarse grids. Directly 
around the shocks finer sub-grid6 are constructed on which we solve the discrete equations after 
adding an artificial viscosity term. We use the grid adaptation techniques developed by Berger [1] 
and Zhu [25]. 

For schemes of even order, Gustafsson and Olsson [5] derived scalar artificial viscosities for con- 
servation laws that support one— or two— point stationary and moving shocks. This theory does not 
directly apply to the Euler equations, but following a similar approach we can construct a scalar 
viscosity that supports near one— point shocks. The viscosity is determined completely by the flow 
variables on either side of the shock, and gives equally good results for the SHOEC-scheme. 

The positioning of the sub-grids is aided by a detection algorithm based on a multi-scale wavelet 
analysis of the pressure grid function which locates potential shocks and spurious oscillations. It also 
supplies the information needed to compute the artificial viscosity terms. The detection algorithm is 
derived from a noise- detection algorithm developed by Mallat and coworkers [13, 14] in signal analysis. 

Most existing shock detection algorithms search for maxima in the first or zero-crossings of the 
second derivative of grid functions. The derivatives are evaluated numerically at the scale of the grid. 
At this scale it is hard to distinguish between maxima that belong to spurious oscillations and those 
that correspond to shocks. Besides, it is difficult to extract information about the shock states when 
the discrete data are distorted by oscillations. If instead the grid function is analyzed on both the 
scale of the grid and a sequence of larger scales, these problems are resolved; the behavior at the 
larger scales is determined by phenomena such as shocks or sharp gradients, while local oscillations 
only influence the smaller scales. The number of available scales is limited because we work in a 
discrete and finite domain. The question is how to construct a reliable algorithm to determine the 
local behavior of the grid function based on this limited information. The wavelet theory supplies the 
answers. In a multi— scale wavelet analysis a function is convoluted with a family of wavelets, each 
varying in position and representing a different scale or frequency. For a special class of families a 
direct correlation exists between the local regularity of the function and the behavior of small, easy— 
to-determine, sets of wavelet coefficients. Moreover, they allow the formulation of a fast wavelet 
transform which leads to an efficient detection algorithm. The underlying wavelet theory is presented 
and the choice of the wavelet family is motivated. We thereby focus on issues that are relevant in the 
context of computational fluid dynamics. 

The sharp shock resolution and the computational efficiency of the above described SHOEC-scheme 
with local grid adaptation are illustrated by numerical experiments for the one-dimensional Euler 
equations. Both shock detection and artificial viscosity are generalizable to multi-dimensions. We 
emphasize that the resulting scheme is a shock-capturing scheme as we do not require any a-priori 
knowledge of the shock location. 

The paper is organized as follows. We start by presenting the Euler equations and relevant consti- 
tutive relations in section 2. Section 3 gives the derivation of the SHOEC-scheme. In section 4 the 
matrix and scalar viscosity are constructed. We give relevant wavelet theory, devise the shock detec- 
tion algorithm and discuss specific implementation issues related to the Euler equations in section 5. 
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In the final two sections we present results and conclusions. 


2. The Euler equations and relevant relations 

The Euler equations in conservation form express the conservation of mass (p), momentum (m) 
and energy ( E ). They are given by 

( 2 . 1 ) u t + f(u) x = u t + f u u r = 0 , 

with state vector u and flux vector / of the form 




p 


pv 

( 2 . 2 ) 

u = 

pv 

, / = 

pv 2 +p 



. E . 


. v(E+p ) . 


The variables v and p are the velocity and pressure of the gets. The latter is related to u through the 
equation of state for a polytropic gas 

(2.3) p = ( 7-1 ){E-^pv 2 ), 

where 7 = 1.4 is the ratio of specific heats c p /c v . The conservation of mass equation, p t — —(pv) T = 
— m x , is often referred to as the continuity equation. The flux vector f(u) is a homogeneous function 
of order one, i.e., f(0u) = 0f(u ). The non-symmetric Jacobian matrix f u is equal to 


(2.4) 


/« = 


0 1 0 
-( 3 ~ 7 )u 2 /2 (3 — 7 ) 1 ; 7 — 1 

(7 - 1) v 3 - *yvEfp 'yE/p — 3(7 — 1 ) v 2 /2 


Its eigenvalues are A] = (v — c), A 2 = u, and A 3 = (v + c), where c = y/ip/p is the speed of sound. 

Viscosity and thermal conduction are neglected in the formulation of the Euler equations. Because 
of the nonlinear and hyperbolic character of ( 2 . 1 ) shock waves can be formed at which these effects 
are of great importance. Additional jump conditions (the so called Rankine-Hugoniot conditions) 
are needed that relate the states Ui and u r on the left and right sides of the shock wave. The 
Rankine-Hugoniot equations for a shock moving with a velocity s are 


(2.5) fi-fr = s(u,-u r ). 

The shock states U\ and u r satisfy the Lax entropy conditions 

A i (u / ) > s > A*(u r ), Ai.^u,) < 5 < A t+ i(« r ), k =1,3. 

In order to simplify the analysis of the equations when there are finite jumps in u and /, Roe 
introduced the linearization 


(2.6) f,-f r = A ,r (Ut-Ur), 

where A lr — f u {u) is the Jacobian matrix in the Roe-variables it, computed such that (2.6) is exact 
[20]. It follows from (2.6) and (2.5) that AjjT = s. The vector u is determined by 

P = \/PiPr , 


and the weighted averages 

y/Pr Vr + y/pi 

y/Pl + yfPr ’ y/Pt + y/fr 

where H is the enthalpy of the system defined by H = (E + p)/p. 


v = 


jj _ y/p^Jb + y7>7 Hj_ 
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The entropy S given by 

S = c„ log (pp~ y ) + constant. 

is used in section 3 to formulate entropy functions of (2.1). In the remainder of this paper we will 
use a simplified scaled expression for the entropy of the form 

( 2 - 7 ) S = log(pp -7 ). 

Smooth flows are isentropic, which means that 


( 2 . 8 ) 


ds os as n 

Dt dt +u dx~ 0. 


where D/ Dt denotes the material derivative. In other words, the entropy is constant along streamlines 
for isentropic flows. Viscous effects and thermal conduction cause the entropy to increase over a shock 
and so equation (2.8) is valid outside shock regions only. 


3. Deriving the SHOEC-scheme 
We consider the system of conservation laws 

(3* 1 ) u t + f x = 0, uJeIR d , xeixoyXi), t> 0, 

u(x,0) = <t>(x)- 

At the boundaries x = x 0 and x = x x of the domain we prescribe data for the in-going characteristics 
in the general form 

( 3 - 2 ) “i(xi,t) = Mi), * = 0 , 1 , 

where u>j are the in-going characteristic variables. 

We employ the energy method to ensure stability. The key to establishing an energy estimate lies 
in a special splitting of the flux derivative vector f z , which we refer to as the canonical splitting , given 
by equation (3.6). To derive it we first write f x as 

(3-3) /* = (/- F). + F x = (f- F) x + F u u x , 

where F = F(u) satisfies Euler’s inhomogeneous differential equation 
(3-4) F u u = -F + f, 

F(0) = f(0), 

and has the form 

(3.5) F(u) = [' f{Ou)d0. 

Jo 

Equations (3.3) and (3.4) now give the desired equation 

(3-6) f x = ( F u u) t + F X = ( F u u) x + F u u x . 

The energy method applied to equation (3.1) yields 

^||u|| 2 = —2 («,/.) = —2 (u, (F u u) t ) — 2 (u,F u u x ). 

The second equality follows from (3.6). If F u is symmetric we obtain 

^|M| 2 = —2uF u u |*J. 

To arrive at a true energy estimate we need to control the boundary term (uF u u), for which we use 
the boundary conditions (3.2). However, these boundary conditions are related to / u instead of F u . 
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In order for the characterstr boundary condition to work we must require that f u and F u have the 
same inertia (i.e., the same number of positive and negative eigenvalues). 

Equation ( 3 . 5 ) shows that F u is symmetric provided f u is symmetric. For the Euler equations this 
is not true however, and thus (2.1) must be symmetrized before the energy estimate can be derived. 

3 . 1 . Symmetrization using entropy functions. Symmetrization of systems ( 3 . 1 ) is possible 
through the use of entropy functions. 

Definition 3.1. A scalar function r] is an entropy function for a system of conservation laws (2.1) if 

• 7] satisfies 

( 3 . 7 ) Tj u f u ~ q u , 

where q(u) is a scalar function called the entropy flux , and rj u = . . . , Tj Ud ] denotes the 

gradient of t/. 

• rj is a convex function of u or, equivalently, T) uu > 0. 

Multiplication of ( 2 . 1 ) by r} u shows that smooth solutions of this system also satisfy the scalar 
entropy equation 

(3.8) T] u u t + Tj u f u u x = 7] t + q x - 0. 

Mock [ 16 ] proved that if an entropy function r) exists for a system of equations ( 2 . 1 ) then the so 
called entropy variables w, given by w T = 7) u , symmetrize the system. This means that u = u(w) and 
f(w) = f(u(w)) have symmetric Jacobians with respect to w. The transformed system is 

(3.9) ut + f(u) x = u t + f(w(u)) x = u w w t + f w w x = 0. 

Harten [6] derived a family of entropy functions using the equation for conservation of entropy ( 2 . 8 ). 
From this equation it follows that 

(3.10) ph(S) t + pvh(S) x = 0, 

for all differentiable functions h(S). If we multiply the continuity equation by h(S) and add equation 

( 3 . 10 ) we obtain 

[ph(S)] t + [pvh(S)] x = 0. 

This is precisely of the form ( 3 . 8 ) with q = ph(S) and q = pvh(S ) . 

The entropy variables w determined by w T = r) u are 

(3.11) w = \e -\ ^—( h/h. - (7 + 1)) -pv p] . 

P L - 1 J 

The state vector u as a function of w reads 

( 3 . 12 ) u = - — r- W3 -to 2 ini - h(h/h - (1 + 1 ))] . 

( 7-1 )h L 1 

The requirement that r] be convex is equivalent to the requirement that w u be positive definite because 
w u = T] uu . Then, the transformation w = w(u) is well-defined as w u is nonsingular. 
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3.2. Energy estimates for the continuous and semi-discrete equations. To obtain the con- 
tinuous energy estimate we apply the canonical splitting to the flux vector derivative f T and the u t 
term. The latter splitting is possible because u w is symmetric. Thus, we define the vector quantities 
U(w), F(w ) € IR d that satisfy the equations 




U w w = -U + u, 

and 



F w w = -F + /. 


As in (3.5) they satisfy 




(3.13) 

U( w) 

= f u(9w ) d9 , 

Jo 

and similarly 

(3.14) 

F(w) 

II 

Si 

¥ 



The above equations show that symmetry of u w and f w implies symmetry of F w and U w . 
The canonical splitting of f x and u, in (3.9) yields 

( U w w) t + U w w t + (F w w) x + F w w x = 0. 

We scalar multiply by w to get 


(3-15) (w T U w w) t + (w t F w w) x = 0, 

and then integrate with respect to x which gives 

(3.16) = - w t F w w |*‘, where 

(3-17) IMItf = / w T U w wdx. 

Jx 0 


Equation (3.13) shows that if u w is positive definite U w is also positive definite, so the norm (3.17) is 
well defined. 

Henceforth, we assume that u(w) and /(in) are both homogeneous functions in w: 


(3.18) 


u(9w) = 0 P u(w) 
f(9w) = 0t>f(w) 


9, (3 em. 


They satisfy Eulers differential equation, i.e., f^w = /3w and u w w = 0u. The homogeneity of / allows 
us to easily relate the right hand side of (3.16) to f w (inertia condition). The homogeneity of u(w) is 
necessary in the derivation of the semi— discrete estimate. For the Euler equations we can construct 
entropy variables that satisfy these homogeneity constraints. 

Equations (3.14) and (3.18) give 


(3.19) 


F(w) = f 1 0 p d9 f(w). 
Jo 


The integral exists if and only if /? > -1 l , in which case we have 


(3.20) 


F(w) = 


1 

l 8 + 1 


f(w), (3 > -1. 


‘Note that a homogeneous function f{w) satisfies Eulers differential equation f w w = (3f. If /? < 0 and 
if / is to be defined for all w in a neighborhood of 0, then / = 0. However, in our applications at least 
one of the elements Wj of the vector w will be positive. Thus we can find a non zero solution of the form 


/(«0 = «£*(*},..., 


W j-1 Wj 


— , . . . , ^ ) for some function g 
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Thus, we can write the canonical splitting of f x as 


(3.21) f x 
Similarly we obtain 

(3.22) 
and 

(3.23) 


(3 + 1 + (3 + J wWt " ft + l ** + (3 f wWx ' 


u ( w ) = 7-T“( w )’ Z 3 > - 1 * 


/3+1 


u ‘ = £7T U ‘ + ^TT U ^‘* 


Although integral (3.19) only exists for /? > —1, the splittings (3.21) and (3.23) are valid for all 
(3^—1 since = / r and u w w t — u t always hold. 

Equation (3.20) leads to 


Furthermore 

IMIfir = jf w T U w wdx = — j-j- jf w T u w wdx = 
The fohowing estimate ca~ tow be shown [18] 


where Ai(tn(x,, <)) are the absolute eigenvalues corresponding to the in-going characteristic variables 
at the boundary in question. As u w and Aj both depend on w this is not an energy estimate in the 
usual sense. But as they are both positive definite we will refer to the inequality as a generalized 
energy estimate . 

When discretizing the equations we apply the canonical splitting only to the spatial variables to 
retain a system of ordinary differential equations (ODEs). Using (3.21) the resulting semi-discrete 
equations are 


(3.25) u, = - (jAj) Df - (j^) /. D», 0*1. 


Here we use bold notation to distinguish discrete variables from continuous variables. Grid vectors 
like u are given by u T = (uq , u[ , . . . , tt^), where tt, 6 IR d , » = 0, . . . , n. The difference operator D 
is defined as 


m J 

(3.26) {Du)j = - j = 0,1,... 

a t=( 

The matrix f w denotes the block diagonal matrix whose blocks are the Jacobians 

It can be shown [19] that u will satisfy the analytic boundary conditions if the initial and boundary 
data satisfy suitable compatibility conditions. In much the same way as in the continuous case, we 
can derive a generalized energy estimate for the semi-discrete equations (3.25) for difference operators 
that satisfy a summation-by-parts principle (B.l) [9, 23]. The second and fourth order difference 
operators that are used in the computations are given in appendix B. The second order operator is 
second order in the interior and first order on the boundary, and the fourth order operator is sixth 
order in the interior and third order on the boundary. 
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3.3. The entropy variables for the Euler equation. The homogeneity requirements of f(w) 
and u(w) restrict the choice of entropy function. Equation (3.12) shows that tz(u;) is homogeneous 
provided the terms h/h and (7 + 1) in the expression for U\ = p have the same dimension. Thus, 
h/h = k, where K is a constant, and consequently 

(3.27) h{S) = Ke KS = K{pp^) K , K, k £ 0. 

This exponential family leads to an homogeneous function u{w) that satisfies u(0w) = u(w) 

which follows from (3.12). Functions h(S ) of the form (3.27) were also discussed by Harten [6]. 
Following his notation we choose k = l/(a+7), a G 1R. It follows that (3 = (a+7)/(l-7). Below it 
is shown that f(w) has the same homogeneity constant. Thus, for a = 1 — 27 the flux vector f(w) is 
homogeneous of order one just as f(u). 

Now, (3.11) and (3.12) yield 

(3.28) w = — [u 3 +- — \p - u 2 u x ] T , 

P 7-1 

(3.29) u = [ «?3 -w 2 w x - - — \p * ] T , 

P 7-1 

where 

(3.30) p* = X e 5/( “ +7) = x (pp-^ a ^\ 

with x = K The variable p* satisfies an equation similar to the equation of state (2.3) 

(3.31) = 

a v 2 w 3 / 

The pressure p can be expressed in terms of p* and w cis 

(3.32) p= i((p*)“^) V(1 - 7) . 

A 

The flux vector 

(3.33) f(w) = ^ [ -w 2 — + p* - — (w! + - — yp*) ] . 

P L w 3 w 3 7 ~ 1 J 

The homogeneity of f(w) follows from equations (3. 32), (3.31) and (3.33). The matrices w u and f w 
are found in Appendix A. The convexity condition on the entropy function given by r) = ph(S ) leads 
to the conditions \ < 0, and a < —7 or a > 0. This excludes the unacceptable case (3 = —1. We 
choose x = “1- 

The equations (3.25) are solved in time by the TVD Runge-Kutta scheme formulated in [22]. For 
the system of ODEs u t = X(ti(tn)), it computes the solution ti n+1 at the (n+l)st time step according 


to 

«i 

= U n 


+ 

AtL(u n ), 

(3.34) 

«2 

= 3/4 u" 

+ 1/4 tii 

+ 

1/4 At 



= 1/3 u n 

T 2/3 U 2 

+ 

2/3Ati(u2), 


with time-step A t and intermediate variables tii and 
Here, L(u(w)) = - (^) Df - f w Dw. 

The variable p* is the most computationally intensive term because of the exponentials in p and 
p. In the discrete equations (3.34) it appears in m, and its reciprocal in f i . For isentropic flow the 
computations simplify; p* is constant and is cancelled. We can use this simplification on the coarse 
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grid, outside the shock regions, where the flow is smooth. We note that for several values of a the 
computation of the exponential terms is fast. For a — 1 — 27 , for example, p* = — p" 2 5 p -3 5 . 

3,4. Non-conservative form. It is well known that if a scheme is in conservative form, the com- 
puted shock speeds will be correct (see, e.g., [10]). This condition is sufficient, but not necessary. 
Tadmor [24] showed how to obtain second-order schemes for initial-boundary value problems that are 
entropy stable and in conservative form. It is however unknown if such schemes exist for higher orders 
of accuracy. The SHOEC scheme (3.25) is not written in conservative form. However, in experiments 
the shock speed has been found to be correct. Below, we show the error in shock location found for 
a single moving shock after long time integration (number of time steps is (9(10 3 )) as a function of 
the step size. It is plotted against the step size. Computations are done in double precision. In all 
experiments, the error was of the order hf 10. Clearly, we can draw the conclusion that the computed 
shock location is correct. 



Figure 1. Error in shock location versus grid size h . 


4. Scalar Artificial Viscosity 


Artificial viscosity is added to the discrete scheme on the refined grids to control the spurious 
oscillations generated around shocks. We derive a scalar viscosity for conservative difference schemes 
of even order that gives approximate one-point stationary and moving shocks. It is determined by 
the flow variables and gives comparable results for the SHOEC scheme. 

4*1. The shock profile. Following the approach in [5] we consider finite difference solutions of (2.1) 
with the Riemann initial condition 


(4.1) 



x < 0 
x > 0. 


The shock states U\ and u r are connected by a k-shock. For simplicity, we furthermore assume that 
U\ is the downstream (supersonic) shock state and that the flow is to the right, so that A k — v — c. 
Applying the coordinate transformation 

y — x — st, 

T = ty 
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(2.1) is changed to the stationary problem 

(4-2) (f-su), = 0, 

= {::>?; 

since u T = 0 in the (y, r)-coordinates. Adding artificial viscosity we discretize (4.2) as 

(4-3) Doifj-suj) = hD + ED_Uj. 

or, for constant E, 

(4-4) f j+i ~ fj-i ~ s ( u j+i ~ u j-i) = 2 E(uj + i — 2uj -f tij-i). 

E is a positive definite matrix. Scalar viscosities correspond to E = el, e > 0, with I the identity 
matrix. The difference operators are defined by 

D ° Uj = 2 ^( u J+ 1_tt i- 1 )’ D + U i = ^(ttj+i-Uj ), D.Uj = 

where h = x J+1 — x : is the uniform mesh size. 

The discrete system is subject to the boundary conditions 

lim Ui = ui, lim Ui = u T . 

i -* -oo * j— +oo 3 

4.2. Supporting approximate one— point shocks. We are interested in matrices E that support 
viscous shock profiles consisting of the states u h u r and an intermediate state u m . The state u m is 
needed in practice to allow for the representation of shocks that are not positioned in the center of a 
grid cell 2 . 

We obtain the following equations at the points x = -h, x = 0, and x = h: 

(4.3) /m ft ~ s { u m ~ u l) — 2i?(u m — ttj), 

(4.6) fr - fi - s(u r - u,) = 2 E(u r - 2 u m + «j), 

(4.7) fr ~ fm ^(n r — u m ) — 2 E{u m — ti r ). 

At all other grid points, the equations are trivially zero. The left hand side of (4.6) vanishes because 
of the RH conditions which gives u m = (u r + u,)/2. As / r = /, + s(u t - u r ), equations (4.5) and (4.7) 
are in fact identical. For the Euler equations however, equation (4.5) can not be satisfied exactly for 
scalar viscosities E = el (see also [7]). Consider the equation corresponding to /j = pv in (4.5): 

(4-8) p,v, - p m v m = (2f + s)(p, - p m ). 

We have p m = (pi + p r )j 2 and p m v m = (p t v t + p r v r )/ 2. Consequently, 

(4-9) piv, - p r v r = (2e + s)(p, - p r ). 

The RH conditions require p\V\ — p r v T — s(pi — p r ), and so 2 c(pi — p r ) = 0, which cannot be fulfilled 
for e > 0. Also, for scalar E equations (4.5) and (4.7) lead to 

(4-10) f m — fi = (s + 2e)(u m -ui), 

(4-11) fr - fm = (s- 2e) (u r - u m ). 

Equation (4.10) represents a A:-shock between states u t and u m moving with speed (s + 2e), while 
(4.11) represents a &-shock between states u m and u r moving with speed (s - 2c). The interaction 
between two such shocks generally results in the formation of a Jb-shock between u, and u r as well as 

2 (4.4) does not support zero-point shocks for E > 0. The left hand side of (4.4) vanishes because of the RH 

conditions. This leads to E(ui — u r ) = 0, which can only be satisfied for uj = u r . 


10 



a contact discontinuity and a rarefaction wave [2] , whence the above described conflict for the Euler 
equations. But, assuming (4.10) and (4.11) are satisfied, a scalar viscosity can be derived that leads 
to approximate one-point shocks. Using Roe-linearization equation (4.10) becomes 

A. m, {u m - u,) = (s + 2e) ( u m - u,), 

where A ml is the Roe matrix determined by the states u\ and u m . Thus, ( u m — Uj ) is an eigenvector 
of A mi with eigenvalue (s + 2c) = A™', or 

c = 1/2(A J"'-s) 

to the left of the shock. Similarly, e = l/2(s - A£ m ) w l/2(\™ 1 — s) to the right of the shock. This 
choice of scalar artificial viscosity corresponds to that used by Jameson in [7]. Since A* > s, A* = s 
and the jump — u m ] is less than the jump [u\ — tt r ], it can be expected that A*" > s also. In 
practical applications we set e = l/2|(Aj” — s) |. We note that IA™* — s\ = min, |A™* — s\. Similar 
analyses for other flow situations lead to the same choice 

(4.12) e = 1/2 min |A™ P - «s|, 

where u p is the upstream state. For scalar equations, the equations (4.5), (4.6) and (4.7) are satisfied 
exactly for the above artificial viscosity. We note that also for systems (4.5) can be satisfied exactly 
if the matrix viscosity E m = ( A lm — si)/ 2 is used. However, E m is not symmetric. So, even if its 
eigenvalues are all positive, E m is not necessarily positive definite. 

The scalar artificial viscosity is often chosen proportional to max,- |A, — s\, Johansson and Kreiss 
[8] show that this viscosity indeed suppresses spurious oscillations, but leads to excessive smearing. 
They also recommend using min, | A* - s |, but evaluate the eigenvalue at u t instead of U/ m , which in 
general leads to a larger viscosity. 

We show the effectiveness of the artificial viscosity (4.12) for stationary shocks in figures 2a and 
2b. They depict the pressure after long time integration of 

(4.13) u t = -D 0 f + hD+ED-u 

using the RK method (3.34) and the standard conservative scheme. The computations are done for 
n = 100, D = Z? 0 , a CFL-number of 0.5 and varying Mach number without local grid adaptation. 
The numerical solutions are indicated by dotted lines and points in the shock region are highlighted. 
The solutions have negligible small amplitude oscillations. 

SHOEC scheme We repeat the computations for the SHOEC scheme with artificial viscosity 
E 

(4.14) u t = - (-j-“j) D 0 f- (j-Ej) f w D 0 w + hD + ED_w. 

The results are shown in figures 3a and 3b, and are similar to those obtained with (4.13). This is 
expected a s f w D 0 w « D 0 f and so (4.14) and (4.13) should be close, at least for weak shocks. We 
remark that if a Roe-linearization fi — f r = A lr (ui — u r ) can be found for the symmetrized Euler 
equations, the corresponding matrix viscosity E m = \A lr — si |/2 would be symmetric. 

Assuming that the shock is located in the interior of the domain, E can be localized to a neighbor- 
hood of that shock; outside this region E would be identically zero. This is possible with a switch as 
defined in (4.21). Hence, 

(4.15) (in, D+ED_w) h — -(D_w, ED_w) h < 0. 

Then, the addition of the term D+ED_w to an entropy conserving scheme will ensure an entropy 
inequality, i.e., the modified scheme is entropy stable. In particular, any scalar viscosity E — el will 
imply (4.15) and thus entropy stability. 
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(b) 

Figure 2. Standard central scheme, scalar viscosity, t/j na j = 2, s = 0; (a) M = 1.1, 
no viscosity (•), scalar viscosity (b) M = 2,scalarviscosity 

Moving shocks Above, the equations are discretized on a grid that moves along with the shock, 
which is not the case in reality. Results for scalar equations are still satisfactory [5]. But, the damping 
by the viscosities is insufficient for the Euler equations as illustrated in figure 4a. This solution is 
computed with 

(4.16) E = 1/2| A lm — si | = l/2(Q|A m ' - sI\Q~ l ), 

where A m ' is the diagonal matrix containing A|"'. For the scalar viscosity e = 1/2 min^ lA™' - s\ an 
identical solution is obtained. We found experimentally that the strong oscillations are primarily 
generated in the A:-th characteristic variable downstream of the shock. This is in agreement with 
the linear analysis by Johansson and Kreiss in [8]. They show that oscillations can appear in all 
characteristic variables upstream of the shock, but only in the k— th characteristic variable downstream 
of the shock, and that the amplitudes of the latter are generally much larger than those of the former. 
The matrix viscosity (4.16) does not introduce much damping on the k— th characteristic variable 
as |A£*' - s\ is small. If we increase the damping by replacing this term with in A m ' in (4.16) 

( I - s | > | AJJ 1 ' - s | ), the spurious oscillations indeed disappear. The result for this adapted 
matrix viscosity is shown in figure 4b. 

We therefore choose the scalar viscosity equal to 

(4-17) e=i|A?'-s| 

for moving shocks. 

4,3. High order of accuracy. We consider the centered discretizations 
(4.18) Q 2r = R 2r D 0 , with 

R*r = 2(-i YPAh 2 D + D-)\ 

v=0 
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Figure 3. SHOEC-scheme, scalar viscosity, tfi na i = 2,5 = 0 ; (a) M = 1.1 ; (b) 

M = 2. 

of order 2 r. The parameters /?„ are determined by 

0o = 1 ? 

0u — t j ^ Pv—h v — 1,2, ...,r 1. 

4/y + 2 

Appropriate approximations of the transformed equation (4.2) with artificial viscosity are 
(4.19) R 2r D 0 (/ - su)j = h R 2r D+ E D-Uj, 

where E and the boundary conditions are as before. In [5] it is shown that (4.19) holds if and 
only if (4.3) holds, and the results produced are the same. To take advantage of the high-order 
approximation property, we must implement the viscosity locally around the discontinuity. We use 


(4.20) R 2r D 0 (f ~ su)j = hR 2r D+ rj. 1/2 ED~ u h 

where the switch r is defined by ([7]) 


(4.21) 


f |A + tif - A_Uf| \ 
VIA+tijI + + 6 / ’ 


6 « 1 . 


Hence, rj = 1 at points where A +Uj and A_ttj are of opposite sign, i.e., at points where the solution 
is oscillatory. In smooth flows r, << 1. We define r^ 1/2 = (r,--! + t 5 )/ 2. By positioning the switch 
and artificial viscosity as in (4.20) the viscosity terms are in conservative form. 


5. Detection with wavelets 

We want to detect shocks and spurious oscillations to aid the local grid adaptation procedure. The 
detection algorithm is based on a wavelet analysis of the pressure grid function. In the first part 
of this section we give a brief introduction to discrete wavelet families (5.1), motivate the choice of 
wavelet family and explain how the local behavior of a function can be determined from its wavelet 
transforms (5.2). The second part deals with implementation issues. The Fast Wavelet Transform 
(FWT) algorithm, introduced in 5.3, makes the detection algorithm computationally attractive. In 
5.4 we discuss the implications of discrete grid functions, finite domains and the presence of spurious 
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Figure 4. Standard conservative scheme, moving shock, initially located at x = 0, 
t/inai = .3 ,Af = 1.5; (a) matrix viscosity (4.16), initial condition (•); (b) adapted 
matrix viscosity (o), adapted scalar viscosity (4.17) (+). 

oscillations. At the end of the section we give a general outline of the detection algorithm. For 
a detailed description of the algorithm, including pseudo— code, we refer to [4]. We recommend 
Daubechies book [3] for a thorough discussion of wavelet theory. The wavelet families used in the 
detection algorithm, and the underlying mathematical theory were developed by Mallat, Hwang and 
Zhong in [12, 13, 14]. 

5*1 • Wavelets and discrete wavelet families. Roughly speaking a wavelet xl>(x) is a well localized 
function that ‘waves’ above and below the x-axis such that / ip(x)dx = 0. 3 Dilations by a E IR+ 
and translations by b £ JR of the mother wavelet ip(x) generate the wavelet family 

(5-2) Mx) = 

The parameter b is the center of the wavelet and a represents its scale. A mother wavelet ip(x) 
and two members of the corresponding wavelet family are depicted in figure 5a. For a > 1 the wavelet 
Tpa.b stretches and can be viewed as representing a lower frequency, while for a < 1 the wavelet narrows 
and represents a higher frequency. Thus, the concepts of scale and frequency are closely related. 

The continuous wavelet transform of a function / e L 2 is given by 

( 5 - 3 ) (fi'/’a.i) = ^(“V~ ) dx > *€lR + ,b€lR 

The variables {/, ip a b ) are referred to as the wavelet coefficients, and give information on the frequency 
content of / near the points b. 


formally, a function rp(x) is said to be a wavelet if and only if its Fourier transform ip(ui) satisfies 


f ! 

Jo W J-oo M 


< +oo. 


This admissibility condition requires the wavelet to have sufficient decay in the frequency domain and that 
V>(0) = 0, which in turn leads to f tp(x)dx = 0. 
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Figure 5. (a) Different members of the wavelet family; (b) Example of discrete func- 
tion /, and the corresponding continuous function f(x). 

For all practical purposes we are, of course, only interested in discrete values for a and b. The 
dilation parameter a is discretized as a m = 2 m , m G Z, where m represents the level. The correspond- 
ing families are referred to as dyadic wavelet families and facilitate the formulation of Fast Wavelet 
Transform (FWT) algorithms (Sec. 5.3) 

At each level m the wavelet coefficients are computed at the centers 6 m ,„ = nAb,n € Z, where A b 
may depend on m. This gives 

X — n./\h 

(5.4) tfrm.n = V> 2 = 2~ m ^( — ), me 2. 

To ensure that the wavelet coefficients contain enough information to characterize /, the family 
consisting of the wavelets V>m,n must constitute a frame [3], i.e., there exist 0 < A < B < oo such 
that, for all / E £ 2 , 

A ll/ll* < E l( /.*»,» >1* < B ll/ll 2, 

m,n€ Z 

The case when A = B = 1 corresponds to m ,„ } being a basis. The choice of the mother wavelet 
V’(i) and the discretization of b are essentially only restricted by the frame and the admissibility 
conditions. The mother wavelet is normally chosen such that it is well concentrated in both the 
spatial and the frequency domain. 

The detection algorithm is based on a wavelet analysis of the pressure grid function given by the 
discrete function /,, with i = 0, . . . , N- 1. Consequently, m and n are bounded; the smallest scale 
that can be considered is determined by the grid size h, the largest scale is determined by the size of 
the domain, and the maximum number of points at which the wavelet coefficients can be evaluated is 
equal to N. Throughout this section, /, is assumed to be the discrete representation of the function 
f(x),x e [x 0 , x jv— i], To facilitate the discussion of the wavelet theory we consider functions f(x) 
that can contain discontinuities at x = x J+i / 2 , represented in the discrete domain by a jump between 
fj and fj+i, and local oscillations at x — Xt when ( fk-i — fk) * ( fk ~ fk- 1 ) < 0. Otherwise f(x) 
is assumed to be slowly varying. An example is given in figure 5b. In section 5.4 we discuss grid 
functions that do not adhere to these simplifications, such as smeared shocks, grid functions with 
spurious oscillations around shocks, and grid functions that vary sharply in the vicinity of a shock. 

5.2. Selection of the discrete wavelet family. The wavelet family is selected based on two 
criteria: 

(1) There must exist a relation between the behavior of the wavelet coefficients and the local 
regularity of f(x). 

(2) The ensuing algorithm must be efficient, i.e., the number of wavelet coefficients needed to 
determine the local regularity of f(x) with sufficient accuracy should be small. 



For very special choices of ip(x') and with A6 — 2 m , the wavelet family VVn.n constitutes an or- 
thonormal basis. Many popular wavelet families, e.g. the Haar, Daubechies or Meyer wavelets, are 
members of this class. The orthogonal families satisfy the first criterion [3], but fail to meet the 
second in this application. To understand this, we observe that shocks and local oscillations may 
travel through the flow domain. The detection results should not change because of their movement, 
in other words, the wavelet transform should be shift-invariant. Figure 6a shows the centers in part 
of the scale-space (a,x) at which wavelet coefficients (/, Vw») are computed in an orthogonal wavelet 
analysis. We assume that f has a bounded discontinuity at x s = 0, at which wavelet coefficients are 
computed at each scale. If the function is shifted by a distance of, for example, 2 ; , only one wavelet 
coefficient is computed at the new location x 9 = 2 J at the levels shown. Also, the wavelet coefficients 
evaluated at positions close to x s change as i, is shifted. In practise, only a finite number of levels 
can be evaluated, and we expect that the detection results, wlrich are based on this limited informa- 
tion, will change as well. Daubechies shows that indeed a very large number of levels is needed to 
accurately compute a for shifted singularities. Specifically, coefficients are needed at scales close to, 
or smaller than, the shift, which may not be available. So, orthogonal families are not suitable for 
detection purposes in the context of this paper; shift-dependency is unacceptable. 


a 


j+3 • < 

• 

j+2 • • i 

• • 

j+1 • • • • i 

9 • • • • 

j 
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0 
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(a) 


a 



"o 

(b) 


Figure 6. (a) Lattice of wavelet centers for orthogonal families; (b) Cone of influence 
of x 0 . 


Clearly, shift dependency can be avoided if the wavelet coefficients are computed at the same points 
at all scales, in other words, if A6 = constant . This constant, which represents the smallest scale 
possible, is determined by the grid size. Because the grid size h is not necessarily equal to a power or 
inverse power of 2, we normalize the smallest scale h to 1 (corresponding to m = 0). Then, A6 = 1 
in (5.4), which gives 

(5-5) m > 0. 

This normalization is used in the remainder of the paper. 

Having thus abandoned orthogonality, freedom in the design of the mother wavelet is gained. The 
question is if this freedom can be used to construct non - orthogonal families that satisfy the two 
criteria given above. Mallat and Hwang [13] showed that this is possible if the mother wavelet rp(x) 
has compact and symmetric support [-K, K ] for K > 0, and if it is a derivative of a smoothing 
function. Henceforth, all mother wavelets are assumed to be in this class. A smoothing function 6(x) 
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is defined as a differentiable function that satisfies 

f d k 0(x) 

(5.6) 0(x) >0, I 0(x) dx = 1, and lim =0, A: = 1, 2, ... . 

From (5.6) it follows that f if>(x)dx = 0 for xj){x) = d k 0/dx k with k > 0. Provided ^(^) has sufficient 
decay in the frequency domain, it will satisfy the admissibility condition (5.1). 

For such mother wavelets the local regularity of f{x) can be computed from a small and easy-to- 
determine set of wavelet coefficients. The local regularity of a function is measured in terms of its 
Lipschitz exponent a, whose definition is given in appendix C. From this definition it foflows that a 
continuous function has a Lipschitz exponent a > 0, a bounded discontinuity (shock) is Lipschitz 0, 
and a Dirac function (local oscillation) has a Lipschitz exponent a = — 1. For ease of presentation, 
the theorem is presented for continuous wavelet families ^ € JR. First we 

make the following observation. If the support of the mother wavelet (p is equal to [-K,K], ip a%b has 
support [b - Ka,b + Ka ]. Thus, at any scale a, the function value f(x 0 ) will influence the wavelet 
coefficients {/, ^p a ,b) computed at the centers b in the set 

(5.7) S{a,x 0 ;K) = {b : \b-x 0 \<Ka}. 

The region C(x 0 ; K) = U a {a) x S(a,x 0 ; K) is referred to as the cone of influence of x 0 (figure 6b). It is 
intuitive to think that the local regularity of / at x 0 should be determined by the wavelet coefficients 
evaluated at centers in C(x 0 ;K). Mallat and Hwang [13] proved that this is true for the functions 
f(x) considered in this paper. Mallat strengthened the result further in the following theorem: 

Theorem 5.1. Let 0 be a smoothing function and let ip k (x) = d k 0(x)/dx k , k > 0 have compact 
support [-K, K], Let b — B(a ) be a curve in the cone of influence of x 0 , with B(a) — > x Q as a — ► 0 
(figure 6b). The function f is Lipschitz a < k at x 0 , if there exists a constant C > 0 and a scale a 
such that along any B(a) 

(5.8) |(/,^ M }| < Ca Q , Va < a. 

Define a modulus maximum at the scale a as a local maximum of |(/, b. The points b at 
which the modulus maxima are attained are referred to as the maximum points at scale a. If f has 
an isolated discontinuity at x 0 , there exist a scale a and a curve b — Bm (a) such that for all a < a, 
Bm{o) art maximum points . Such curves are referred to as maximum lines. The Lipschitz exponent 
of f at x 0 is determined by (5.8) evaluated along the maximum line converging to x 0 . 

A similar theorem also holds for dyadic families, with a m = 2 m , and b £ IR [15, 14]. However, we 
cannot compute coefficients at scales smaller than 1, nor at positions other than the grid points. Thus, 
we presume that the decay of the available wavelet coefficients at scales larger than 1 characterizes the 
regularity of the function. If this is not the case, the mesh must be refined (decreasing h ). The second 
part of the theorem makes the wavelet detection computationally attractive; the Lipschitz exponent 
a at a point x 0 can be estimated from a limited number of wavelet coefficients. At a sequence of 
levels m, where 1 < m < m mar < log 2 (A r ), we compute the wavelet coefficients at the grid points. 
Next, we find the modulus maxima at each level and locate maximum lines that converge on the scale 
of the grid. Then, we determine a such that 

(5.9) log 2 |(/,V£, n )| = log 2 C + ma, 

as close as possible, in a least squares sense, for the wavelet coefficients along the select maximum 
lines. 

As an example, figure 8b shows a signal with several points of distinct behavior, and the wavelet 
coefficients computed with the detection algorithm at increasing levels + 1,.... If the Lipschitz 
exponent a is positive, the amplitude of the modulus maxima increases with j. The smooth variation 
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points 1 and 4 exhibit this behavior. The singularity at 3 produces modulus maxima that decrease 
as the scale increases, so sup(o) < 0. The amplitudes of the modulus maxima corresponding to point 
2 do not change and we conclude that sup(a) = 0. 

The parameter k limits the maximum Lipschitz exponent that can be determined. As we are 
interested in detecting discontinuities (a = 0) and spurious oscillations (o < 0), we choose k = 1. 
We remark that for k = 1, a direct connection exists between the wavelet detection and detection 
algorithms that search for maxima in derivatives. Define 0 Oi& = a _1 0( £ ^). The function 0 a t is a 
smoothing function also as / 8 a h (x)dx = f 6{x)dx = 1. Furthermore, 


so that 


i’a.b = a 


M ajL 

dx ’ 


(5-10) ^ ) = -«!< /. «... >. 

In other words, the wavelet transform is the first derivative of the function f smoothed at the scale 
a by 0 a ,b- So, detecting modulus maxima corresponds to detecting maxima in derivatives. It fol- 
lows that a discontinuity has one maximum line converging to its location as the derivative of the 
smoothed function will have an absolute maximum at this position. A Dirac function naturally has 
two maximum lines. Furthermore, we note that the wavelet coefficients provide an estimate of the 
jump [/] a at scale a across the shock: 


(5.11) 


l[/]-l 


I [/] a I 


a 


db 


= l< 


which can be used as an extra criterion for shock selection. 


5.3. The Fast Wavelet Transform. For special the wavelet transformations can be per- 
formed through fast discrete filter operations [4, 14]. Given two sequences {</*}* and {h™} k , we 
define G m = {g ™ }, and H m = {h™}, as the sequences constructed by putting (2 m_1 -l) zeros between 
the elements of, respectively, {jr*}* and {h k } k . Also, W m = {(/, Tp m ,n)}n is defined as the sequence of 
wavelet coefficients at level m. It can be shown that certain rf{x) possess sequences and {h k } k 

such that the wavelet coefficients at the levels ni ^ 1 can be computed by the discrete convolutions 


(5.12) 


W m 


X~Gm * <S m _ 1 

A m 

H m * 5 m _] 


m = 1,2, .. . 


where the sequence So contains the grid function values /;. The sequences S m represent smoothed 
versions of / at level m, whereas W m can be seen as representing the details lost in this smoothing 
process. For this reason, the sequence {h(k)} k is referred to as the impulse response of a low pass filter 
and {<7(1:)}* as the impulse response of a high pass filter. The finite resolution of f(x) necessitates the 
scaling by the parameters A m . If the sequences S m were available at all levels m < 0 (for a m — ► 0), the 
scaling parameters would not be needed. We refer to [4] for full details concerning their derivation. 
At each level m, we perform O(n) operations. The maximum number of levels is equal to log 2 (n). 
For detection purposes it is sufficient to compute the wavelet coefficients at a total of four or five 
levels. Therefore, the overall detection algorithm is O(n). In the detection algorithm we use the 
mother wavelet depicted in figure 8b. It is a quadratic, anti-symmetric spline with compact support 
[-1, 1] given by 


(5.13) 


ip(x) 


2(x + l) 2 

-1 < x < 

-4x(l + x) - 2x 2 

-1/2 < x < 

— 4x(l — x) + 2x 2 

0 < x < 

— 2(x - l) 2 

1/2 < x < 


- 1/2 

0 

1/2 

1 
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Its filter coefficients h k and g k , and the scaling parameters A m are given in table 7. Note that 
symmetry is obtained by applying a shift of 1/2 to W m . 


Filter Coefficients 

Scaling parameters 

k 

h k 

<7* 

m 


m 


-1 

0.125 


i 

2 

5 

1.336 

0 

0.375 

-2.0 

2 

1.5 

6 

1.334 

1 

0.375 

2.0 

3 

1.375 

7 

1.333 

2 

0.125 


4 

1.344 

> 7 

1.333 


Figure 7. Filter coefficients and scaling parameters for the spline wavelet 



Figure 8. (a) Example wavelet transformation; (b) The spline wavelet ip 


5.4. Implementation issues. The pressure grid function /, may be distorted by spurious oscilla- 
tions around the shock. They are characterized by wavelength L,L >2h and amplitude r. Figure 9a 
shows the damping factors r///r and h h l t ■, where r // and r // u are the amplitudes of the smoothed 
oscillations after respectively one and two applications of H . Oscillations with a wavelength L — 2h 
are removed immediately, and oscillations with short wavelengths are dampened quickly. Therefore, 
the oscillations have little and negligible effect on the computation of the Lipschitz coefficient a if we 
use levels m > 2 in (5.9). An example is given in figure 9b. It shows a pressure distribution obtained 
from a central finite difference scheme applied to the Euler equations. The oscillations are prominent 
at the smallest scale but quickly disappear. Note that the jump |[/]| across the shock is accurately 
represented by the wavelet coefficients. 

If a shock is smeared, the estimated Lipschitz exponent will be larger than zero. However, as the 
scale increases, the modulus maxima will quickly converge to a constant, equal to the total jump over 
the smeared shock. For one- and two-point shocks, the correct Lipschitz exponent can generally be 
detected from the modulus maxima at the levels m > 3. Otherwise, local oscillations and/or large 
magnitudes of the wavelet coefficients at the larger scales will indicate that there is a steep gradient 
necessitating local grid refinement. We remark that as a result of local grid refinement and an efficient 
artificial viscosity, smearing is generally not an issue at the coarser scales. 
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Figure 9. (a) Damping Factors Filter H, and (b) Wavelet coefficients {/, V'm.n) cor- 
responding to an oscillatory grid function /*. 

If the grid function varies sharply close to a discontinuity, the Lipschitz exponent will also be 
perturbed. However, the grid adaptation procedure will refine the grid locally as the detection 
algorithm will pick up local oscillations. At the finest scales, the variation in the grid function close 
to the shock will not be as strong and the discontinuity can again be detected with sufficient accuracy. 

As the discrete domain is finite, we need to periodize /,. To avoid the introduction of artificial 
shocks at the boundary, the finite discrete domain is extended by reflection according to 

( 5 - 14 ) fi = hn-i, for n < i < 2n-l. 

Although the reflected image is continuous at the boundary, its derivative may not be so. At the 
boundary itself this will not cause any problems because the filter G is anti-symmetric. In the local 
grid adaptation procedure, we ensure that the shocks are located in the interior of the (refined) grids, 
and that the flow varies smoothly near the boundaries. Therefore, the cusp or kink resulting from 
a discontinuous boundary derivative will be very small and will have negligible effect on the wavelet 
coefficients in the domain interior. , 

Approximations to the shock states U| and u T are needed to construct the artificial viscosity. If no 
local spurious oscillations are present, this task is trivial once the shock location is known. Otherwise, 
the data at the scale of the grid is distorted and we approximate /t and f r from the smoothed data 
S m at scales m > 1. 

We summarize the findings concerning the choice of the wavelet family: 

• The family is dyadic for computational purposes. 

• The family is redundant to obtain shift-invariance. 

• The mother wavelet is compact, differentiable and the first derivative of a smoothing function 
so that theorem 5.1 is valid. 

• A Fast Wavelet Transform exists, leading to an O(n) detection algorithm. 

Below, a general outline of the detection algorithm is given with references to the relevant equations, 
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theorems and sections. 


A general outline of the detection algorithm 

• Periodize the grid function /, by reflection. 

Reference 
Eq. (5.14) 

• Compute {/, ipm, n) in the extended domain at 
the scales a m = 2 m , m = l,...,m mor with a 
Fast Wavelet Transform. 

Eqs. (5.12), (5.12), 
Table 7 

• Locate the maxima of | (f,ip m ,n)\ at each 
scale. 


• Determine the converging maximum lines. 

Thm. (5.1) 

• Compute the Lipschitz exponents from the 
evolution of the wavelet coefficients along the 
converging maximum lines. 

Eq. (5.9) 

• Locate shocks, sharp gradients and spurious 
oscillations. 

Sec. 5.2, 5.4 


6. Numerical results 


High order of accuracy 

The initial conditions 


( 6 . 1 ) 


u(x, 0) 


( 1 + 0.5cos(87rx) 
tii(x,0) 

^T+2 u i( x ’°) 


with periodic boundary conditions, correspond to the linear advection of a cosine density wave. The 
velocity and pressure are equal to 1 and remain constant in time. This test problem is chosen to 
illustrate the advantage of difference operators with a high order of accuracy. The flow is smooth and 
therefore solved without local grid adaptation or artificial viscosity. Figure 10 depicts the solutions 
at t = 0.3 obtained for the SHOEC scheme with the sixth and second order difference operators given 
in the appendix 4 , and the first order Godunov scheme for 50 grid points. 

To achieve the same accuracy as the sixth order scheme, a minimum of 250 grid points is needed 
for the second order scheme and at least 2000 grid points for the first order scheme. 

Shock detection, artificial viscosity and high order of accuracy 


4 Note that because of the periodic boundary conditions the global accuracy is sixth and second order 
respectively as well 
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Figure 10. Solid line exact solution, 6th order 2nd order (•), 1st order (-). 
We consider the initial condition 


( 6 . 2 ) 


u(x,0) 


/ 3.857143 \ 

10.141852 ] x < -4 
\ 39.166666 / 

C l + 0.2 sin 5x \ 

0 ) x > -4 

2.5 ) 


A strong shock wave, moving to the right, interacts with a sine wave. of small amplitude. The 
shock speed s = 3.55, resulting in a Mach-number Af « 3 in front of the shock. The difficulty lies 
in resolving the low and high frequency oscillations that are formed behind the shock wave in the 
density. This problem is a good test for resolution, as well as the detection algorithm and artificial 
viscosity. The solution is not known analytically. For comparison, an "exact” solution is computed 
using a formally conservative 4-th order central scheme with n = 2000 (no grid adaptation) and 
t = 1- The conservative form will ensure that the shock location is correct. The density is shown in 
figure 11a. Again, the computed solution is indicated by a dashed line and the points in the shock 
are highlighted. The shock is resolved very sharply which illustrates the effectiveness of the artificial 
viscosity. 

Figures lib and c show the coarse grid solutions for a fourth order SHOEC scheme (sixth order 
in interior and third order at the boundary) obtained for n = 200 and n = 68 respectively. We used 
local grid adaptation and one level of refinement (n is the number of coarse grid points). The shock 
location is determined with the detection algorithm; a priori knowledge of the shock speed is not 
used. The ratio of fine to coarse grid points is 4 for n = 200 and 8 for n = 68 . We used the 4th order 
SHOEC scheme whose difference operator is given in appendix B (6th order in interior, 3rd order on 
the boundary). In figure lid we depict the location of the fine grid relative to the coarse grid as time 
progresses for n = 200. The true and detected shock locations coincide. Clearly, the refined grid is 
confined to the shock region. The variations in fine grid size are due to spurious local oscillations 
away from the shock that are detected and included in the fine grid. The solution for n = 200 is very 
accurate. For n = 68 the fine grid solution is good and on the coarse grid the principal features of the 
flow are clearly visible. The average number of coarse and fine grid points in this last computation 
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is 150. The computation took approximately 8% of the time for n = 2000. The detection algorithm 
and the regridding procedure only accounted for 5% of the total time. 




(a) 


(b) 



Figure 11. Density (a) conservative scheme, n = 2000; (b)+(c) local grid adaptation 
SHOEC scheme, coarse grid (-), fine grid (•), coarse grid point (o), fine grid point (+), 
true shock location (*); (b) n = 200; (c) n = 68; (d) Position of refined grid, n = 200, 
true and computed shock position (*). 


7. Conclusions 

We have presented a high order scheme for the Euler equations that satisfies a semi-discrete 
energy estimate, as well as an efficient strategy for computing sharp shock solutions. The high order 
of accuracy of this Split High-Order Entropy Conserving (SHOEC) scheme leads to a very efficient 
scheme for flows that are smooth in large regions of the flow domain. The construction of the SHOEC 
scheme is based on a symmetrization of the equations, followed by a canonical splitting of the flux 
vector, and difference ope.ators that satisfy a summation by parts principle. The same techniques 
can be applied to other systems of conservation laws that satisfy the cone condition, mentioned in 
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section 3. We are currently developing a SHOEC scheme for mixed hyperbolic-parabolic equations 
in ocean modeling. 

Although the scheme is in a non-conservative form, the shock location obtained in numerical 
experiments is correct and converges with h. This is illustrated in figure 1. In [21], Salas and Iollo 
formulated shock jump conditions for the Euler equations in their primitive, non-conservative, form. 
Jump conditions for the symmetrized split Euler equations are being investigated. Their existence 
may lead to an analytical explanation for the correct shock speed. 

We constructed a cheap and effective scalar artificial viscosity that is used in locally refined grid 
in shock regions. It results in sharp shock resolution. The positioning of the finer sub-grids is aided 
by a detection algorithm based on a multi-scale wavelet analysis of the pressure grid function. The 
algorithm detects shocks and local oscillations efficiently. No a-priori knowledge of the shock location 
is used and the resulting method is therefore shock capturing. The obtained information on the shock 
position is used only in positioning the sub-grids and constructing the artificial viscosity. It could 
prove to be useful in other applications as well. As an example, we mention the high-order methods 
by LeVeque and Shuye [11] that need an accurate shock location. 

The sharp shock resolution and computational efficiency of the above described SHOEC-scheme 
with local grid adaptation are illustrated by the numerical experiments. 

The SHOEC scheme, artificial viscosity and shock detection are generalizable to multi-dimensions. 
In a forthcoming paper we will present numerical results for multidimensional applications in which 
we use the technique of overlapping grids as well as local grid adaptation. 


Appendix A. Symmetrization 

For h(S) = e s ^\ 


where 
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Here, p* and p are related through equations (3.30) and (3.32), R = l/(a + 7) and R = (R- 1)(7~ 1). 
For a = (1 — 27), in is a homogeneous function of degree 1 in u, and 
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The symmetric matrix f(w) w is given by 

Cm Cpv 2 — p 

Cpv 2 -p v(Cpv 2 -3p) 

[ v(f pv 2 - § pv< - p 

where C = (1 — a — 7 )/or. 




v(Zpv 2 -^) 

f pv 4 -^£-(l + 


7-1 


pv z 


^ - (I + )P v2 ^ P ** + Jpv 4 ) J 


Appendix B. Difference Operators 
We use difference operators that satisfy the summation by parts 

(B.l) (u ,Dw) h = u^w n - u^w 0 - (Du,w) h , 

for arbitrary discrete grid vectors u and w , with respect to the weighted scalar product 

n 

(B.2) = h (TijuJwj. 

i,}= 0 

The scalar product can also be written as 


(B.3) 


(u,w) h = h u t Ew. 


To establish the semi-discrete estimate for the Euler equations, E must be diagonal [19]. In [9], 
Kreiss and Scherer showed that there exist diagonal E and difference operators D of accuracy 2 p 
in the interior, 1 < p < 4, such that the summation-by-parts property (B.l) holds. The difference 
operators D and matrices E are computed in [23], 

The simplest example is the second order accurate difference operator 


with 



-1 1 

-0.5 0 0.5 


-0.5 0 0.5 

-1 1 




A sixth order accurate difference 


operator D with third order boundary closure is given by 
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j = n — 5, ... ,n 
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The values of d,-,-, 0 < j < 5,0 < t < 8 are given below. The corresponding S is diagonal, i.e. 


0, i ^ j. Its diagonal elements 

On = 

Oi are 
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0 

^30 

= 
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^31 

= -30637/64308 
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0 

^34 
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72/5359 
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See [23] for a derivation of the coefficients. 


Appendix C. Lipschitz exponent 

Definition C.l (Lipschitz exponent). Let 0 < a < 1. 

(1) A function /(x), defined on a domain T)j C 7R , is Lipschitz ot over an interval 

( x h x r ) C V } ,if and only if there exists a constant C < oo such that for all x,y € (x h x r ) 


!/(*) “ f(v ) I < C \x - y|°. 


(2) A function /(x) is Lipschitz at at a point Xo, if and only if there exists a constant C < oo such 
that 


sup 

*€X>/ 


!/(*)- /(*o)l 

|z - X 0 |° 


< c. 


(3) A function that is not Lipschitz 1 at a point x 0 is said to be singular at that point. 

(4) A distribution is Lipschitz a over an interval (x ; ,x r ), if and only if its primitive is Lipschitz 
at + 1 on (xt,x r ). 

(5) A distribution f(x) has an isolated singularity Lipschitz a at a point z 0 , if and only if f(x) is 
uniformly Lipschitz a over an interval (xi,x r ) with io € (xi,x r ), and /(x) is Lipschitz 1 over 
any subinterval of (x,,x r ) that does not include x 0 . 
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